Division of Engineering 
BROWN UNIVERSITY 
PROVIDENCE, R. I. 



FINITE-ELEMENT FORMULATIONS 
PROBLEMS OF LARGE ELASTIC- 
PLASTIC DEFORMATION 

R. M. McMEEKING AND J. R. RICE 



04 - 2841 ©^ 
Onclas 

G3/^2 4 33 90 J 

Technical Report 
NASA NGL 40-002-080/15 to the 
National Aeronautics and Space Administration 


'32 p tLi ? c IIC MFOBHSTIOH (Bpobb Oniv.) 

CSCL 20k 


May 1974 



Finite-Element Formulations for Problems of Large Elastic -Plastic Deformation 


by R. M. McMeeking and J. R. Rice 

Division of Engineering, Brown University, Providence, Rhode Island 

May 1974. 


Abstract 

An Eulerian finite element formulation is presented for problems of large 
elastic-plastic flow. The method is based on Hill’s variational principle for 
incremental deformations, and is ideally suited to isotropically hardening 
Prandtl-Reuss materials. Further, the formulation is given in a manner which- 
allows any conventional finite element program, for "small strain" elastic- 
plastic analysis, to be simply and rigorously adapted to problems involving 
arbitrary amounts of deformation and arbitrary levels of stress in comparison 
to plastic deformation moduli. The method is applied to a necking bifurcation 
analysis of a bar in plane-strain tension. 

The paper closes with a unified general formulation of finite element equa- 
tions, both Lagrangian and Eulerian, for large deformations, with arbitrary 
choice of the conjugate stress and strain measures. Further, a discussion is 
given of other proposed formulations for elastic-plastic finite element analysis 
at large strain, and the inadequacies of some of these are commented upon. 

* 
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Introduction 

Many elastic-plastic problems of interest involve large deformations and 
the present paper is concerned with finite element formulations for this class 
of problem. Such formulations are also relevant to small deformation elastic- 
plastic problems. Indeed the conventional "small strain" formulation in stress 
analysis can be inadequate for small strain problems not only when increments of 
rotation greatly exceed those of strain (buckling of slender members) but also 
when representative stress levels have attained a magnitude comparable to that 
of the plastic hardening modulus. The latter can occur after only minute deforma- 
tions for lightly hardening metals. These situations relate more to the shape of 
the member under consideration and to intrinsic material properties, rather than 
to the smallness or largeness of strain per se, but are best approached in a rig- 
orous finite deformation context. 

Several workers have proposed and utilized finite element schemes for problems 
of finite deformation, and these include diverse formulations of both the Lagrangian 
(material) and Eulerian (spatial) type. The next section gives a brief review of 
these and includes some comments on ambiguities in certain of the approaches to 
elastic-plastic behavior. Thereafter we proceed directly to an Eulerian formula- 
tion which is ideally suited to the large deformation of Prandtl-Reuss materials. 
This is given in a manner which allows the straightforward adaptation of any exist- 
ing small strain program to finite deformation analysis. It also serves as a basis 
for analysis of the possible inadequacies of the conventional small strain formula- 
tion. The program is applied to a numerical study of the initiation of plastic 
bifurcation in a plane strain tensile bar. 

Finally, in the concluding section of the paper we have presented a general 
formulation of the finite deformation problem in finite- element terms. By appro- 
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priate specialization this can be made to coincide with various separate formula- 
tions, both Lagrangian and Eulerian, and provides a suitable basis for comparison 
of the different approaches as well as for a fuller discussion of ambiguities in 
other elastic-plastic formulations as noted above. 

Review of Finite Elements for Finite Deformation 

The description Lagrangian is attached to large deformation finite element 
programs which use a mesh of elements representing some fixed reference state for 
strain. A scheme of this type was proposed by Hibbitt, Marcal and Rice [1], who 
derive their finite element rate equilibrium equations from the principle of 
virtual work for large deformation. They identify four stiffness terms, which 
are named small strain stiffness, initial load stiffness, initial strain stiffness 
and initial stress stiffness. In elastic-plastic analysis, all of these must be 
calculated for each increment of deformation and the last three have a complicated 
form. Needleman [2] has a similar Lagrangian scheme, but he derives his equations 
from a variational principle due to Hill [3]. A third Lagrangian formulation was 
presented by Felippa and Sharifi [4], who intend to place no limitation on the size 
of an increment of deformation. Higher order terms are thus introduced into their 
stiffness. However, these terms are not significant in a tangent modulus approach to 
elastic-plastic or non-linear elastic analysis, since the increment of deformation 
must be small in any case to ensure that the constitutive rate moduli do not change 
significantly from one increment to the next. For the materials mentioned the 
method of [4] coincides with that of [13. 

The Eulerian formulation which is based on a mesh representing the current 
state of deformation has been used by Yaghmai and Popov [5]. In their presentation 
they derive equations from a current configuration version of the variational 




principle used in [4] and solve an elastic problem. However, Sharifi and Popov 

[6] extended the method of [5] to elastic-plastic analysis of infinitesimal 
strains but finite rotations , although they do not seem to consider the possi- 
bility that plastic deformation moduli may be of a size comparable to current 
stress levels. Another Eulerian formulation is due to Gunasekera and Alexander 

[7] , who base their equations for a Prandtl-Reuss material on the principle of 
virtual work. The finite element equilibrium equations they obtain are correct 

if their stress o is interpreted as the 2nd Piola-Kirchhoff stress, but a state- 
ment to this effect is not made. If their o is so identified, the constitutive 
law used does not then seem to constitute an adequate form of the Prandtl-Reuss 
equations. There also seems to be an ambiguity in work on an Eulerian formula- 
tion for Prandtl-Reuss materials by Argyris and Chan [8], who do not make a clear 
definition of their stress and strain terms. Two interpretations appear possible, 
but neither constitutes an adequate finite element method. This point is discussed 
in the final section of the paper, along with further comments on [63 and [73. 

Osias [9] has an Eulerian scheme, which admits non-symmetric constitutive laws 
through a Galerkin method. He derives the same rate equilibrium equations as ob- 
tained here, but his constitutive law leads to a non-symmetric stiffness. 


An Eulerian Finite Element Formulation for Finite Deformation 


A concise formulation of rate equilibrium at arbitrary amounts of deformation 
is given by the following form of the virtual work equation, which has been cited 


by Hill [3]. 
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where all integration extents are in the reference configuration and X is the 

position vector of a material point in that reference state, t is the non-sym- 

metric nominal stress, which is defined so that the force vector f per unit 

reference area of a surface having the reference unit normal n° is given by 

f. ~ n?t.. , and b is the body force per unit reference volume. Rates are 
3 i i] 

indicated by the superposed dot and 6v is an arbitrary virtual velocity varia- 
tion which disappears where velocity rates are prescribed i.e. on S° - S° , 
where S° is the reference surface on which tractions are prescribed. Now, by 
the device of choosing a reference state that is instantaneously coincident with 

the current state, t may be simply related to a spin-invariant stress rate, the 

* 


Jaumann or co- rotational rate of Kirchhoff stress t , more suited to use in 
constitutive relations. The relationship is [3] 


hi * T ij " “kjhi ' “iAj + a ik T j ,1c 


( 2 ) 


where o is the Cauchy stress, v^ ^ is 3v^./9x^ where x is the position 
vector of a material point in the current state and D is the rate of deformation 
tensor such that 
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The Kirchhoff and Cauchy stresses are related through 


(3) 
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(4) 


where J is the ratio of volume in the reference state to volume in the current 
state (so that J = 1 instantaneously). 

Under these circumstances (1) becomes [3] 
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f . 6v.dS + b.Sv.dV 
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• • 

where all integration extents are in the current configuration, f and b are 
still nominal force intensity rates » but with respect to the current areas. 

The principle (5) may be used to formulate an Eulerian finite element method. 
For this purpose we define {$>} , which is the vector of rates of nodal degrees 
of freedom and {v),(e} and {D} , which are respectively the vectors of compo- 
nents of the velocity, the infinitesimal strain rate and the deformation rate 
all at a point in an element. In a conventional small strain program the position 
and element- type dependent matrices [N] and £B] relate {v} and {e} to the 
rates of degrees of freedom of the undeformed mesh so that 

{v} = DOU> , t e } = [B]{$,} . (6) 


However in an Eulerian program, the same matrices [N] and [B] , as used in a 
conventional small strain program for elements of a given kind, relate {v} and 
{D> to the rates of degrees of freedom of a mesh representing the current 

geometry, i.e. 

{v> = [N]{5»} , {D> * [B]{i>} . (7) 


The matrices are related by 




( 8 ) 


where [N^J , [B„3 are rows of DO , [B] defined so that 


v^ * . 


The finite element rate equilibrium equations obtained from (5) have two 
stiffness terms. The first stiffness arises from T..6D..dV , and the other 

j Y J-J 

is an initial stress stiffness [k ] arising from the remainder of the left-hand side 

s 
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of (5). Note however that there is no initial strain stiffness, [k 3 is obtained 

s 

from the relationship 


3U> 

s 


26D, . c « » D, * ) dV , 
ki i] kj 5 


or 


* } * / ¥ < 4 v k,i °ij Vj - 


(9) 


where { 6 ^ } is the variation in nodal rates of the current mesh that generates 
6v , 5D . Hence the coirplete finite element equations are 


f 


CB] [C][B]dV + [k s ]l {+} = {P> , 


L 


( 10 ) 


where {P> « I CN] T {b}dV + f [N] T {f} dS . 


[C] is the rate independent incremental constitutive matrix appearing in the 
form 


(t > = [C3 {D} 


( 11 ) 


and {x } is a vector of components of t . When the prescribed force rates 
are such that {b} and (f) on S T cannot be fully identified, {P} divides 
into a prescribed part and a geometric part, linear in {t|i} . The last part 
must be added as another stiffness term in the manner described by Hibbitt, 

Mar cal and Rice [13. 


The form (5) given above and thus this finite element method is of general 
validity. Although the constitutive law of interest 


T ij = ^ijki D k* 


* 


( 12 ) 


or equivalently {t } = [C] {D} is materially objective to superposed spins, 

the tensor £ in its most general form may be dependent on the total deformation 
from some prior state including total rotations. However (5) and the associated 
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finite element method is of greatest convenience for problems in strain rate 
independent isotropically hardening elastic-plastic materials with rate con- 
stitutive laws written in terms of t and D and dependent only on material 
parameters and the state of stress. This class of materials includes those with 
different loading and unloading paths. For these materials the finite element 
rate equilibrium equations for the current mesh are independent of the total 
deformations. Furthermore the first stiffness in this case is formed identically 
in terms of constitutive rate moduli as in small strain programs. The second 
stiffness has a simple universal form in terms of a . Indeed, the nature of 
these stiffnesses allows a convenient adaption of small strain programs with an 
appropriate constitutive law. 


Constitutive Laws 


The constitutive law used to formulate any finite element equations at finite 
deformation may lead to a non- symmetric stiffness. The most general case leading 
to a symmetric stiffness is a law derived from Hill's [33 homogeneous quadratic 
rate potential 4> , so that 


and 


*ij = a(av]/3x.) = BOv^/ax^wav^)^) 

♦ 1 k »y«i • 


(13) 


Such laws are symmetric and rate independent and have been phrased in terms of 
conjugate variables of stress and strain. All elastic constitutive laws with 
work-potentials and elastic-plastic laws with a normality rule expressed in terms 
of conjugate variables of stress and strain take this form [3,10,11]. The 
elastic-plastic laws have zero order dependence on strain rate, since the response 
at yield is governed by the direction of loading. However, when both the elastic 
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and the elastic-plastic moduli are applied to a strain rate tangential to the yield 
surface, the same stress rate results. It is then obvious that at yield (13) is 
still uniquely satisfied and * has continuous first variation with respect to 
strain rate. For materials for which $ exists, (1) may be rewritten as 

S { f v o *<^> dV ° - | v o Vi dV ° - L Vi dS °} 1 0 (14 > 

and a current configuration version of (14) is derived from (5) so that 

‘ jj v M ? ,dV - 1 - v m v m mv -j v V i dV - j s "i v i dS } = 0 > 

T (15) 


where 


u * Kj l - D) 


and (16) 

* 3U _ 3 2 U y> 

T ij “ 3D. . " 3D.. 3D,. k£ ijki ki ’ 

J 13 13 ki J 

Hill [12] has studied the transformation between constitutive matrices for dif- 
ferent stress rate and strain rate measures and has also shown [ 3 ] that the 
existence of <t> for a material implies the existence of U , and of similar 
quadratic potentials when constitutive laws are phrased in terms of conjugate 
stress and strain measures; the converse is so also. 


Prandtl-Reuss Equations 

The classical Prandtl-Reuss equations for isotropically hardening materials, 

as presented e.g. by Hill ([13], pp, 15 to 39), are not restricted to small-strain, 

although it is necessary to specify for them a suitable materially objective stress 

ft 

rate when principal deformation axes are rotated. The Jaumann rate 0 of Cauchy 
stress coincides with a in the absence of rotation and hence may be chosen for 
their generalization. The equations are then 
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_ 1+v * v . * 

D ij " e °ij " e 5 ij kk 


for elastic loading or any unloading and 


1+v * 


0 » » * 6 • * 0 \ \ + 

E 13 E 13 kk 


9o! .a 0. # 
13 ki kl 


for loading at yield, where 


°ij " 3 d ij °kk * 


s 2 a lj a lj * 


v is Poisson's ratio, E is Young's modulus, h is the slope of the Cauchy 
stress - lograithmic plastic strain curve for a simple tension test and 6^ is 
the Kronecker delta. The equation (18) may be inverted to show that for plastic 


loading 


* 

a. . * 

13 


1+v C6 ik 6 j* l-2v S ij S kl- 2 -2 ( 4^ } ] * 


There is no rate potential for this law, and when Osias used it as a generalization 
of the Prandtl-Reuss equations he obtained a non- symmetric stiffness. Hibbitt et 
al. used the same generalization, and Budiansky [143 has proposed laws which can be 
shown to be equivalent to the equations (18,19) given above [15], [16]. Budiansky's 
laws , which use Green strain and which are referred to curvilinear co-ordinates , 
were used by Needleman [2] and Chen [17]. A generalization of the Prandtl-Reuss 
equations, which does admit the potential U and does lead to symmetric stiffness 




E ft r + ^ r Jt ^ _ 1 T\ 

1+v t6 ik 3! + l-2v 13 ki „-2,2 L . E , J D k* 




for plastic loading and 


T ij = 1+v C6 ik 6 jA + l-2v 6 ij 6 ka ] D k4 
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for elastic loading or any unloading. This law is referred to a current config- 
uration and is equivalent to (17), (18) and (19) to within terms of stress 
divided by elastic modulus since volume change is taken to be purely elastic. 
Needleman and Chen note this equivalency and although they postulate that 
Budiansky's equations govern material response, they use equations equivalent to 
(20) for their calculations. This last form has also been introduced and dis- 
cussed by Hutchinson [163. It should be noted that the elasticity in (20) is 
approximately derivable from a work potential to within terms of stress divided 
by elastic modulus in comparison to unity [18]. As stated previously a rate 
potential must evidently exist for elastic response, and (20) but not (19) is in 
accord with this. The matrix [C] for the law (20) is homogeneous of degree zero 
in {J) , in the usual manner, for elastic-plastic materials at yield. Since 
this [C] is also the constitutive matrix relating stress rate to {e} in a 
small strain program, the stiffness of such a program calculated in current con' 
figuration constitutes one term of the stiffness of the Eulerian finite deforma- 
tion program. To this must be added [k 3 to form the complete stiffness. 

S 


Conventional Small Strain Formulation 


The conventional small strain principle of virtual work may be written as 


a . .ct.dV = f b.vtdV + [ f.vtdS 

J v ^ ij J v ii 'S x 1 


( 21 ) 


where a is the stress field in equilibrium with surface tractions £ and body 
forces b , is the virtual rate of the infinitesimal strain tensor and is 
compatible with the virtual velocity field v^ , S is the surface of the body 
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and V is the volume of the body. All these functions are associated with the 
undeformed state and the difference between deformed and undeformed volume is 
neglected. In that case, since the virtual field v* may be taken as time inde- 
pendent and also taken as some arbitrary variation 6v , vanishing wherever v 
is prescribed on S (i.e. on S-S^) , 


o. ,<$e..dV 

13 ij 


b.6v.dV + 

J v i i 


L 


f^6v^dS . 


( 22 ) 


Suppose that the constitutive rate law has the form 


°ij * *^ijk£ £ ki 


(23) 


where X = X (a) is symmetric in interchange of ij,ki. X takes different 
forms for plastic loading versus unloading, but o varies continuously with c . 
Then (22) is equivalent to Cl9] 


61 [v] = 0 , where 


(24) 


u ? ] ■ L - L Vi" - 1, 


Vi* 


and this variational principle may serve as the basis for a finite element formu- 
lation. Alternatively and equivalently , (22) could serve as well and would be 
needed for non-syrranetric • Using the previously defined terms {tji} , {v} , 

{ c } , [N] and [B] of (6) for small strain analysis, and rendering I of (24) 
stationary for the discrete class of fields generated by {ip} , leads to finite 
element equations in the usual way. 

In the introduction certain inadequacies of conventional small strain anal- 
ysis were mentioned. These may be illustrated by considering equation (15) and 
assuming now that the current deformed state is indistinguishable from the unde- 
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formed state. In this case all stress measures are identical, but all stress 

rate measures are not. These differ by terms of order o times 3v/dx . It 

* * 

seems plausible to think of o in the small strain formulation as being t , 

since this may be assumed to satisfy the postulated Prandtl-Reuss law. Under 
these circumstances (15) becomes 



1 

2 


o. .(e. .)e. .dV 

i] i] 



c^.(2e 


ik e kj 


7wr )dV - / v Vi* 


3v, 3v. 
3X 



0 


(25) 


The additional terms proportional to o that are neglected in a small-strain 
formulation [compare (24) with (25)3 are well-known to be important for slender 
structural members, in which rates of rotation can greatly exceed rates of strain- 
ing. 


However, even when buckling is not a possibility, it is evident that the 

neglected terms of the small-strain formulation are of the same inportance as 
* 2 * * 

those of order he arising from the ce part of (25), wherever h has a 
magnitude comparable to that of current stress levels. In this circumstance, 
predictions of, say, the slope of the overall load-displacement curve in the 
plastic range cannot be considered accurate. This does not invalidate the solu- 
tions from small-strain programs , but means only that one should not consider 
that the actual stress-strain relation in the hardening range has been precisely 
modelled. Indeed, the approximation is of the same kind as when the ideally 
plastic idealization is adopted, with the yield stress taken as some representa- 
tive flow stress for the amount of deformation at hand. Of course, h is often 
far greater than current stress levels and then, in the absence of buckling tend- 
encies, small-strain programs are unobjectionable. 
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Some Numerical Results 

The Eulerian elastic-plastic finite deformation scheme was incorporated in 
the manner described into an existing small-strain program, that is a modified 
version of the relevant part of the MARC program developed by Professor P. V. 
Marcal, and in use on the IBM 360-67 computer at Brown University. This program 
models elastic-plastic Prandtl-Reuss materials in a manner described by Tracey 
[20]. An iterative procedure attempts convergence to a best solution during an 
elastic-plastic increment, incremental stiffness being calculated by the partial 
stiffness approach of Marcal and King [21], subsequently modified by Rice and 
Tracey [22], 

A first examp le illustrates the treatment of finite rotations by this scheme, 
and uses 4 elements forming a square, which was stressed to represent a state of 
elastic plane strain. Of the 10 degrees of freedom, 8 were constrained to cause 
a rotation of ir/4 radians in first SO and then 25 increments. The results 
correspond closely to the exact solution as given by a simple Mohr circle trans- 
formation, and are shown in the table. 

A solution for the initiation of tensile necking of an elastic-plastic bar 
in plane strain with a constant hardening modulus was also obtained. Osias [9] 
studied a similar problem. The bar had a ratio of length to width of 3 and the 
center l/6tb of the length was thinned by 0.5% so the analysis could be formu- 
lated as a standard deformation problem rather than an eigenvalue problem for 
abrupt bifurcation. One quarter of the bar was analyzed using finite elements 
and this section is marked OABC in fig. 1 of the specimen with the thinning 
greatly exaggerated. There were 432 elements in the mesh and 241 nodes, as shown 
in fig. 2. Uniform strain elements in the crossed array shown were chosen because 
Nagtegaal, Parks, and Rice [23] have recently shown this array to be exceptionally 
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free of artificial constraints for incompressible or near incompressible deforma- 
tion. The material modelled had a ratio of Young's modulus to equivalent Cauchy 

2 

yield stress of 7.5 * 10 , Poisson's ratio was 0.3 and the ratio of hardening 

modulus to yield stress was constant at 1.25. The results are shown by the full 
line in fig. 3, where the amplitude 6 of the neck (i.e. , the width reduction at 
the center section), normalized by the initial length L , is plotted against 
engineering strain y = AL/L . A Fourier analysis of the neck showed that the 
lateral bar surfaces deformed most dominantly in cosine shape that is predicted 
from the rigid plastic analysis of Cowper and Onat [24], In fig. 3, point a 
is maximum load, point b is the point at which elastic unloading first occurs 
and point c corresponds to the specimen illustrated in fig. 4 with the boundary 
between the elastic elements and the plastic elements marked zz » Those ele- 
ments , which have unloaded elastically, are above the line zz . 

It is known [25] that abrupt plastic bifurcation in an elastic-plastic or 
rigid-plastic body must occur with the overall deformation increasing at a def- 
inite initial rate with respect to the bifurcation amplitude. This rate just 
causes neutral local stress alteration at a single point or locus of points, with 
subsequent local unloading toward an elastic or rigid state. The dashed line in 
fig. 3 is drawn at a slope in accord with that of the initial strain vs. bifurca- 
tion amplitude slope predicted from the Cowper-Onat solution. The line is located 
so as to cross the vertical axis at the critical strain predicted from their re- 
sults (phrased in terms of a critical ratio of hardening rate to true stress, 
dependent on the current length to width ratio for the bar) . It is seen that the 
large growth of the initial imperfection as predicted from the finite element 
solution is in close accord with their results. Cowper and Onat found that neu- 
tral loading occurs at a point corresponding to point C in fig. 1, and suggested 
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that unloading would spread from there towards the side AB of the bar and down 
to the center line OA . This occurred in the numerical analysis of the elastic- 
plastic bar, and it can be seen from fig. 3 that once elastic unloading occurred, 
this analysis began to deviate from dashed line initial slope for the rigid- 
plastic bar. 

The numerical results do not, however, show the inception of unloading at 
strains very near the amplitude at which large growth of the imperfection sets 
in, as would be expected from the general result noted earlier for abrupt plastic 
bifurcations. Also while it is clear that concentration of deformation in the 
neck and elastic unloading at the end of the bar are both occurring in the numer- 
ical analysis, neither are occurring to a great extent. Another aspect of the 
analysis is that the elements in the neck are elongated in one direction, which 
has reduced their accuracy. Any further extension of the bar would only aggravate 
this problem, and thus it would not be profitable to continue the analysis. The 
elongation of these elements may seem to reduce the value of this Eulerian formu- 
lation as opposed to a Lagrangian formulation, where the elements retain their 
shape. However, the same effect is embedded in the stiffness of the Lagrangian 
elements and they have no advantage in this respect. 

The latter inadequacy at large deformation in the neck is clearly traceable 
to element size; and smaller elements may also cause the inception of neutral 
loading to take place earlier. For example, Osias [9] used a finer mesh and, 
while no comparison of his results with the Ccwper-Onat solution was made, he 
did find the inception of unloading at a lower strain and a faster spread of the 
unloaded zone down the bar axis. 
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Finite Element Formulations for Arbitrary Deformation Magnitudes 

Let e denote any material strain tensor. This is defined on a coordinate 
system X in the adopted reference state such that the principal directions of 
e coincide with the fibers of principal stretch from that state, and so that 
its principal values e^ , e ^ , e^ are related to the corresponding stretch 
ratios Aj , A^ > ^y a ,nono ‘ ton i cai ly increasing function g( A ) : 

e s g(A) , with g(l) a 0 and g 1 (1) = 1 . 


Any strain measure so defined accords with the infinitesimal strain e at small 

* 

displacement gradients, and e = D initially. The most commonly used of this 
type is Green strain 

e = y (A - 1) , e ij = | (u i,j + u j,i + u k,i u k,j ) * (26> 

where u = x - X is the displacement from the reference state and, now 
u. . = 3U./3X. . We shall also have occasion to use logarithmic strain 

e^ a log A (27) 


although this cannot be simply expressed in terms of 3u/3X . However, by series 

G 

expansion of an arbitrary measure e in tern® of e , one obtains 

e a e G + (m-l)(e G ) 2 +...., where m = ~ [g n (l) - 1] ; (28) 

G L 

the parameter m = 1 when e = e , and m = 0 when e = e . From this an 
arbitrary strain is given to second order in 3u/3X by 


5 . . = e. . t — u. .u. . + (m-l)e., 

ij 2 kjiTc,] lk k] 


. , where e.. - 4 (u. .+u. .). (29) 

13 * t ,3 ] j 1 - 


Following Hill [11], we define a corresponding family of symmetric conjugate 
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stress measures s , such that s..6e.. is the stress working per unit volume 

x] 13 

of reference state for arbitrary virtual deformations 6e . For example the 

PK G 

second Piola-Kirchhoff stress s is the conjugate to e , where 



3X. 3X. 

J K v 


3x 

J = det(^) . 


(30) 


As suitable starting points for incremental finite element formulations, 
either Lagrangian or Eulerian, one may proceed directly from the rate equilibrium 
form (1) and transform t to the rate measure most convenient for the constitu- 
tive description as done earlier or, alternatively, proceed in the manner of [1] 

Q 

from a virtual work statement in conjugate deformation variables (taken as e 
PK 

and S in Clj). The end results are the same, and to illustrate this we now 
take the latter approach, starting from 

f »e) T {s)dV° * f {Su} T tb)dV° ♦ f («u) T (f}dS° , (31) 

V J V° J S° 

where <5u and 6e are associated virtual variations, the forces are nominal, 

and s , e are any conjugate measures . Matrices [N] and [B] generate the 

6 quantities from variations in the nodal degrees of freedom vector: 


Uu} * [N] {6*} , (fie) * [B] {6vp} . (32) 

Moreover, [B] is now a function of 3u/3X , which is different for each differ- 
ent strain measure, although all coincide with a common [B°] , the same as that 
in (6,7,8), when the current and reference states coincide. Thus the finite 
element equilibrium equations are 

f [B] T {s}dV° = {P> , (P) = [ [N ] T {b}dV°+ f ^ [N] T {f)dS° , (33) 

J v ° Jv° 

and the incremental form is 
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< CB3 T {s> + [B3 T {s})dV° = {P} 


where [B] will be a linear measure dependent function of {J;} , and where 
{P} is given by (10), except for being evaluated in V and on S , at least 
in the normal case as considered here for which the interpolation function [N3 
has no dependence on {i|»} . 

A more precise identification of [B3 and [B] is obtained by first noting 

from (32) that a given column of [B] can be identified as 3{e}/3i{i where 

a a 

is some member of {i|/} . Thus if we regard (e) as a function of 3u/3X , and 
recall the row vectors [N^] composing [N] , 

[b: . |ili_ CO . (35) 

This allows [B3 to be written as 

[5] . <») T CV!t [ V,j • (36) 

The finite element equilibrium equations that follow from (34) are thus 

(O c ] ♦ Ck s 3){5)} = {P} , (37) 

where the incremental stiffness arising from the constitutive matrix [C] in 

(s) = [C] (e) (38) 


■ J v o 


CB] T [c3CB)dv° . { D.,3* f£i- BJ m- CN k ) it dV< 


and that arising from the existing stress state (i.e. , the initial stress 
stiffness) is 


fv TOT ‘ s, ^^ dv ° • 


Ck s ] ^ 
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Both of these are measure dependent, although in a compensatory way, because 
[C] and the derivatives of (e) are measure dependent. 


G 

When the adopted measures are Green strain {e } and second Piola-Kirchhoff 
stress {s PK > , one finds from (26,32) or (35) that a given row » 
corresponding to the variation <$e^ , is 




(41) 


and further the initial stress stiffness is 


tv 






(42) 


This seems to be the most suitable basis for a Lagrangian method, since few 

other strain measures are readily computed at finite strain. It is important, 

however, that the constitutive matrix [C] be properly identified. For example, 

ft ft 

if the Prandtl-Reuss relations are taken to hold in terms of (t } or (o } and 
{D} , then the matrix [C°] in 

{t*} = [C°]{D} or {a*} = CC°]{D) 


must be properly transformed to CC] , as for example in [1] . 


By contrast, Eulerian finite element formulations are readily made for any 
choice of measure by taking the reference state to coincide instantaneously with 
the current state. Then (s) = (a) and (e> = (D) momentarily and [B] has the 
small deformation form 


[ V * T [ V,j + T c V.i 


( 43 ) 


when expressed in terms of the current geometry of the element. By setting 
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(s) = {c} in (40) and evaluating the mixed derivatives of {e} from (29) when 
3u/9x = 0 , the initial stress stiffness is 




■I, 




+ 2(m-l)[B, .1 a. .[B. ])dv 
Ki 13 3* 


(44) 


Although the stress measures (s) and {a} coincide, their rates do not 
and once more special attention must be paid to the matrix [C] in the stiffness 

[k ] = I [B] T [C][B]dV . (45) 

C ' V 


The importance of proper identification of [C] in both the Lagrangian and 
the Eulerian approaches is clearly indicated when the reference state is taken 
to coincide instantaneously with the current state. Then as Hill [11] has shown 


E « ■ t i 3 - ■‘■’uAj * D ikV * ( “ 6) 

where m is still the measure dependent coefficient of (28). Since m = 0 

when the adopted measure is logarithmic strain, this means that the Jaumann 

ft 

( co- rotational) rate t of Kirchhoff stress coincides instantaneously with the 
rate of the stress measure conjugate to logarithmic strain. 

o 

Thus, suppose that is the constitutive tensor relating t to D . 

Then the constitutive tensor for any other choice of measure when the current and 
reference states coincide is given from (46): 

T ij = * ijkAl i "' ,Ues E ii = <47) 

In the special case m = 0 , [k 3 and [k 3 of (44) and (45) are then the stiff* 

s c 

nesses given earlier in equations (9,10). 

Equation (47) without m specifically identified shows the error that would 
be entailed in writing the Prandtl-Reuss equations in the form (20), but using 
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• 9? 

some arbitrary measure s rather than t as the stress rate. This step be- 
comes permissible as an approximation only when the smallest components of , 

having the order h , greatly exceed components of a , the magnitude of which 

*PK 

may be represented by o . In particular the use of s (m=l) would imply 

in uniaxial tension or compression a plastic hardening rate of h + 2o when 

phrased in terms of true stress vs. logarithmic strain. If o is of order h 

this is in conflict with the approximately symmetric tensile and compressive 

response of initially isotropic materials (when phrased in true stress and 

logarithmic strain) which the equations are intended to represent [13]. This 

is pivotal to an assessment of refs. [6,7,8]. In the first two the constitutive 

law (20) is formulated with s taken as the stress rate, although in [7] a 

definition or description of the adopted measure of stress is not given. However, 

the terms [B] and [k ] in [7] are clearly computed for Green strain, which 

s 

requires that the stress rate used must be the rate of Piola-Kirchhoff stress for 
validity of the equilibrium equations. The problem concerning the constitutive 
representation is not severe for [6], its intent being only to treat small strains. 
In that case h often far exceeds the absolute value of o and the imprecision 
as to stress rate is unimportant. However the adopted constitutive representa- 
tion is not suitable in [7], which is intended for large plastic deformations, 
when the absolute value of a will generally be greater than h . 

Argyris and Chan [8] do not define what measures of either stress or strain 
are intended in their work. They cannot be inferred unambiguously since no ex- 
plicit formulae are given for their term analogous to [B] . Indeed it is not 
clear that their stress measure is a member of the conjugate family, but for what- 
ever measure it is , the Prandtl-Reuss equations are written with its ordinary 
time derivative as the stress rate. If their form for [B] is intended as for 
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Green strain, then their procedure is identical to that of [7] and all comments 
concerning [7] also apply to [8]. On the other hand if their stress is meant to 
be Cauchy or Kirchhoff stress, then the use of its ordinary rather than co-rota- 
tional time rate in the Prandtl-Reuss equations creates a constitutive law that 
is not spin invariant. 
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Table 1. 

Finite Rotation by 45° in the x-y plane for a stressed plane-strain specimen. 
Comparison of 50 and 25 increment finite element solutions with exact Mohr-circle 
trans formation . 



°xx 

a 

yy 

0 

xy 

Before Rotation 

1.000 

0. 

0. 

Mohr Circle 

.500 

.500 

.500 

50 increments 

.500 

.500 

.513 

25 increments 

.501 

.499 

.526 


Figures 


Figure 1 

Figure 2 
Figure 3 

Figure 4 


Specimen for necking analysis, with imperfection shown greatly 
exaggerated. The actual thinning is 0.5% of full width. 


Finite element mesh for necking analysis. 


Engineering strain y versus amplitude of the neck 6 normalized 
by the undeformed length L . 


Deformed mesh at engineering strain 0.69. Elements above line zz 
have unloaded elastically. 
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Figure 1 
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Figure 4 


